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Abstract 



The paper develops a solver based on a conforming finite element method for a 3D-1D 
coupled incompressible flow problem. New coupling conditions are introduced to ensure a 
suitable bound for the cumulative energy of the model. We study the stability and accuracy 
of the discretization method, and the performance of some state-of-the-art linear algebraic 
solvers for such flow configurations. Motivated by the simulation of the flow over inferior vena 
cava (IVC) filter, we consider the coupling of a ID fiuid model and a 3D fluid model posed in a 
domain with anisotropic inclusions. The relevance of our approach to realistic cardiovascular 
simulations is demonstrated by computing a blood flow over a model IVC filter. 
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1 Introduction 

Coupling a 3D fiuid fiow model and a system of hyperbolic equations posed on a ID graph is a 
well established approach for numerical simulations of blood flows in a system of vessels [32] . Such 
a geometric multiscale strategy is particularly efficient, when the attention to local flow details 
and the qualitative assessment of global flow statistics are both important. The relevance to 
cardiovascular simulations and challenging mathematical problems of coupling parabolic 3D and 
hyperbolic ID equations put 3D-1D flow problems in the focus of intensive research. Thus, the 
coupling of a 3D fluid/structure interaction problem with a reduced ID model merged to outflow 
boundary, which acts as an absorbing device, was studied in |15j . The coupling of a 3D fluid problem 
with multiple downstream ID models in the context of a finite element method was considered in 
|38) . In [33], a system of a 3D fluid/structure interaction problem and a ID flnite element method 
model of the whole arterial tree was implemented to model the carotid artery blood flow; and in [S] 
a unified variational formulation for multidimensional models was introduced. A splitting method, 
extending the pressure-correction scheme to 3D-1D coupled systems, was studied in [26] . 




In most of these studies, the 3D model was a generic fluid-elasticity or rigid fluid model, while 
numerical validations were commonly done for cylindric type 3D domains (with rigid or elastic 
walls) ; several authors considered geometries with bifurcation |38[ 133) or constrained geometries 
(modeling a stenosed artery) [^. More complicated geometries occur in simulations of blood flows, 
if one is interested in modeling the effect of endovascular implants, such as inferior vena cava 
(IVC) filters. In numerical simulations, a part of a vessel with an intravenous filter leads to the 
computational 3D domain with strongly anisotropic inclusions. A downstream flow behind the 
implant may exhibit a complex structure with traveling vortices, swirls, and recirculation regions 
(the latter may occur if plaque is captured by the filter). Moreover, the IVC flow is strongly 
influenced by the contraction of the heart, and both forward (towards the heart) and reverse (from 
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the heart) flows occur within one cardiac cycle. Downstream coupling conditions for such flows 
may be a delicate issue. Thus, the flow over an IVC filter is an interesting and challenging problem 
for a 3D-1D flow numerical solver. 

The coupling conditions of 3D and ID fluid models and their properties were studied by several 
authors. Coupling conditions and algorithms based on subdomain iterations were introduced in 
|15j . and the stability properties of each subproblem were analyzed separately. The first analysis 
of two models together was done in [T5]. In that paper, it was noted that if the Navier-Stokes 
equations are taken in the rotation form and suitably coupled with a ID downstream flow model, 
then one can show a bound for the joint energy of the system. It is, however, well known that 
using a finite element method for the rotation form of the Navier-Stokes equations needs special 
care |22j . and setting appropriate outflow boundary conditions can be an issue. In the present 
paper, we introduce an energy consistent coupling with a ID model for the convection form of the 
Navier-Stokes equations. The joint energy of a coupled 3D-1D model is appropriately balanced 
and dissipates for viscous flows. 

Handling highly anisotropic structures is a well-known challenge in numerical flow simulations 
and analysis. There are only a few computational studies addressing the dynamics of blood flows in 
vessels with implanted filters. Recently, Vassilevski et al. [34| numerically approached the problem 
of intravenous filter optimization using a finite-difference method on octree cartesian meshes to 
resolve the geometry of implants. In that paper, it was also discussed how the effect of an implant 
can be accounted in a ID model through a modification of a vessel wall state equation (see also 
[36l[37] for the development of this method for atherosclerotic blood vessels). In the present paper, 
we take another approach and locally resolve the full 3D model, while keeping the state equation 
unchanged. We report on a finite element method for modeling a 3D-1D coupled fluid problem, 
when the 3D domain has anisotropic inclusions. Naturally, this leads to meshes containing possibly 
anisotropic tetrahedra. We study the performance of the finite element method both by considering 
the accuracy of solutions and by monitoring the convergence of one state-of-the-art linear algebra 
solver for the systems of linear algebraic equations to be solved on every time step of the method. 
We are interested in the ability of the solver to predict such important statistics as the drag force 
experienced by an intravenous implant. 

The remainder of the paper is organized as follows. In section [2| we review 3D and ID fluid 
models and discuss coupling conditions. The stability properties of the coupled model are also 
addressed in section [2] Section [3] presents a time-stepping numerical scheme and an algebraic 
solver. In section|4j we validate the 3D finite element solver and the coupled method by considering 
the benchmark problem of a flow past a 3D cylinder and a problem with an analytical solution. 
The application of the method to simulate a blood flow over a model IVC filter is given in sectionjsj 
Numerical experiments were performed using the Ani3D finite element package |40j , which was used 
to generate tetrahedra subdivisions of 3D domains, to build stiffness matrices, and to implement 
the linear algebra solvers described in section [3) 

2 The 1D-3D coupled model 

This section reviews 3D and ID fluid models and describes the coupling of the models. In this 
study, the 3D model is assumed 'rigid'. 

2.1 The 3D model 

Consider a flow of a viscous incompressible Newtonian fluid in a bounded domain £7 C M^. We 
shall distinguish between the inflow part of the boundary. Fin, the no-slip and no-penetration part 
(rigid walls), Fq, and the outflow part of the boundary, Fout- On the inflow part we assume a given 
velocity profile. The outflow boundary conditions are defined by setting the normal stress tensor 
equal to a given vector function 4>. Thus, the 3D model is the classical Navier-Stokes equations in 
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pressure- velocity variables: 



du 

~dt 



+ (u • V) u - i/Au + Vjj = 



r!x (o,r], 



u|ri„ = Ui, 
du 
9n 



div u 
u|ro 



: 

0, 



(1) 



Here n is the outward normal vector to 917. The system is also supplemented with initial condition 
u = Uq (div Uq = 0) for t = in fi. 

We remark that the notion of 'inflow' and 'outflow' boundary is used here and further in the 
text conventionally, since the inequalities u-n<Ooru-n>0 are not necessarily pointwise 
satisfied on Fin or Foutj respectively. In applications we consider, the mean flux, u • n averaged 
in space and in time, is expected to be negative at Tm and positive at Tout- However, for certain 
t G [0,T], the flux J u • nds may take positive values at Fin (negative at Fout)- 
ri„(out) 

If the solution to ([I]) is sufficiently smooth and the inflow boundary conditions are homogeneous, 
the following energy balance holds: 



^-||u|P 
2dt" " 



I'llVull 



u|-^n - </) • uds = 0. 



Here and in the rest of the paper, |j • || denotes the L^(r2) norm. If one assumes 



|u|2u-nds>0 Vte[0,T], 



(2) 



then solutions to ([!]) satisfy the a priori energy inequality and this opens possibilities for showing 
partial well-posedness results. Even though, the assumption ^ was used for the purpose of analysis 
in the literature on 3D-1D blood flow models (see [121 [S]), it is hard to verify ^ for practical 
flows. Moreover, the inequality ^ no longer holds if reverse flows occur, as, for example, happens 
inlVC [mill]. 



2.2 The ID model 

A one-dimensional model can be derived from the Navier-Stokes equations posed in a long axisym- 
metric elastic pipe by integrating over cross section, making some simplifying assumptions and 
considering integral average quantities as unknowns, see, e.g., [Tl[32]. Let ui{t,x) be the cross sec- 
tion of the pipe normal to x, S(t, x) is the area of ui{t, x) and u{t, x) is the axial velocity. Introduce 
the averaged variables: the mean axial velocity u and the mean pressure: 

u{t,x) = S~^{t,x) / u(t,x)ds, p{t,x) = S^^{t,x) / p(i,x)ds. 



We consider the model given by the following system of equations for unknowns S: 

dS d{Su) 



dx 

du d{u^/2 + p/p) 



(p{t, X, S, u) 



= 'i/'(i, X, 5, u) 



for X e [0,1]. 



dx 



(3) 



For initial conditions, the mean velocity and the cross section area are prescribed, u|t=o — ^Oi 
S\t=o — Sq. Here Pcxt is the external pressure, ip{t,x, S,u) is a function modelling the source or 
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sink of the fluid, as may be required in hemodynamic simulations, if a blood loss happens in a 
vessel. Further, we assume tp — and p^xt — 0, so from now p has the meaning of the difference 
between the fluid pressure and the external pressure. The term ip(t,x, S,u) accounts for external 
forces, such as gravity or friction. Following PQHH], we set 

t/j ^ ^16iyurj{S){Sd^)-\ S ^ S-^S. (4) 

Here v is the viscosity coefficient, d is the pipe diameter, S is the reference area (in the hemodynamic 
applications S is the cross section area of a vessel at rest) and 



2, for S > 1 
S + S-\ forS^<l. 



The last equation in ([s]) relates the pressure to the cross section area. The function / is defined by 
the elasticity model of the pipe walls, cq is the elasticity model parameter. We use the one from 

m- 

f expfS-S^^i - 1) - 1, forS'>S^ 
fiS)^\ L ^ (5) 

[ In(S'S-i), forS'<S'. 

Other algebraic defining relations linking the mean pressure and the cross section area are known 
from the literature, see, e.g., [32]. They are equally well suited for the purpose of this paper. 

In [15], the authors derived the energy equality for the one-dimensional fluid model in Q,S,p 
variables, where Q = Su is the flux. Up to possibly different choices of 'ip{t,x, S,u) and f{S), the 
formulation in [15j is equivalent to (|3]). Written in terms of u, S, and p, the energy equality for ^ 
is (recall that we assumed ip — 0): 



d 



1 



^^^Dit)-pJ S'ip{t,x,S,u)udt ^ - Su{p+ -u ) ^ (6) 



with the energy functional 



£iD{t)^^l SuM^ + pc^ f f /(s)dsdx. 



For f{S) given in ([5|, the second term in the definition of £iij{t) is always positive, making £i£)(t) 
positive for all t > 0. The choice of ^{t, x, S, u) in Q ensures that the second term on the left-hand 
side of (|6| is positive as well. Thus, for the homogenous boundary conditions the energy of the ID 
model dissipates: ^Sioit) < 0. 

System ([s]) is hyperbolic and can be integrated along characteristics. To see this, we write ^ 
in the divergence form: 

dV dF{V) _ 

-? - 
u p 

with V = {S, u}, F = {Su, 1 — }, g = {0, Tp}- Denote by and (1=1,2) the left eigenvectors 

2 p 

and eigenvalues of the Jacobian A = One finds 



. _ , u+i-iycoyS/SeMS/S-l) ifS>S „ 

A,; — < ' , I — i, Z, 



and 

T 



co^{i/s) exp{s/s-i), i-iyVs 



ifS> S 

, i-1,2. 



(co, (-1)^5)^ if S'<S^ 
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Figure 1: The schematic couphng of JIsd and JI^'d™" domains. 



Thus, system ([s]) can be written in the characteristic form: 



dV_ ^dV 

dt * dx 



i = 1,2. 



(7) 



Under physiological conditions in hemodynamics, it holds u < cq, implying that the eigenvalues 
Xi have opposite signs. Therefore, Wi is the incoming characteristic from point x = 1 and W2 is 
the incoming characteristic from point x — 0. Two boundary conditions, one at a; = and one at 
X — 1, are enough to close the system. 



2.3 The coupling of ID and 3D models 

We now consider two domains J^sd and ri^'^^", as shown in Figure [l] In fi^J^™" we pose the ID fluid 
model and in fl^D we pose the full three-dimensional Navier-Stokes equations. The one-dimensional 
domain fJ^'g"" is coupled to the downstream boundary of fisD • 

There are several options to define the coupling conditions of ID and 3D models, as discussed, 
for example, in |15j . One can ask for the continuity of the mean pressure, the mean axial velocity, 
the flux, the normal cross section area, the averaged normal stress, or the entering characteristic. 
In general, the continuity of all these quantities cannot be satisfied simultaneously. A choice has 
to be made. 

One common choice, see, e.g., |331 15] . is to impose the continuity of the normal stress and the 
flux: 



/ du 
\ dn 



pn 



P\x=dn, 



/ u- nds = Su\x=d- 



(8) 
(9) 



This choice is, however, known to be energy inconsistent in the following sense. Assume the 
homogeneous boundary conditions on Fjn and the downstream end of r^^*™". Then the cumulative 
energy balance of the coupled model is 



dt 



i£3D{t) + £iD{t)) + iy\\Vu\\ 



-{p+ ^\u\'^)n] ■ uds + Su{p+f^u^' 
on 2 



(10) 
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with S-ioit) — |||u|p. K^{S) is a positive coefficient defined from Q. Easy to see that for couphng 
conditions (|8|-|9]) the right hand side of ( 10 1 reduces to 



P 
2 



Su'' 



luPfu. 



n)ds 



In general, it is not clear if this quantity is non-positive and thus if the cumulative energy is properly 
dissipated, as it holds for the full 3D Navier-Stokes equations with the Dirichlet homogeneous 
boundary conditions. To circumvent this inconsistency, it was suggested in jl6j to replace condition 
(Is]) by the continuity of so-called total stress: 



dw. , 

V— + [p - 



/- , P-2 



(11) 



9n ' ' 2' 

Condition ( 11 ) together with the continuity of the flux, i.e. condition ([9|, makes the right hand side 
of ( 10 1 equal to zero. Hence, the cumulative energy dissipates. However, setting the total stress 
equal to a constant is not a consistent outflow condition for the simplest Poiseuille flow. Moreover, 
(111 is the natural boundary condition for the rotation form of the Navier-Stokes equations. Using 



it with the common convection or conservation forms leads to non-linear coupling conditions and 
often requires iterative treatment [7]. Although the rotation form is an interesting alternative to 
the standard convection form, it is still not a standard option in the existing software and its use 
requires certain care |22j |23] . 

The condition for the normal stress, as in (|8]), is the natural boundary condition for the com- 
monly used convection and conservation forms of the Navier-Stokes equations. Such a condition 
has been shown to be surprisingly useful as outflow boundary condition 19 . Thus, instead of 
keeping (|9| and changing ([s]) to the total stress condition, we retain (|8| and instead of ^ assume 
the continuity of the linear combination of the fluid flux and the energy flux: 



p{d) j u-nds+^ j |u|2(u-n)ds 



(12) 



One easily checks that the combination of (|8|) and ( 12 1 ensures the right-hand side of ( 10 1 equal 



to zero and thus the correct energy balance and energy inequality are valid as formulated in the 
following theorem. 

Theorem 2.1 Consider the coupled 1D-3D fluid problem given by ([l]), ([s]), Q, and coupling 
conditions (Is]), (12 1. Assume the homogeneous boundary conditions Ujn = and u\x=e = 0. Then 



a sufficiently smooth solution satisfies the following energy decay property: 

£3D{t)+£lDit) + l^ f ||Vu||2dt'+ / / K,{S)uMxdt' ^ £3d{0) + £id{0) (13) 

JO JO JO^i™" 

for any t € [0,T]. If tfj is defined m then K^{S) = 16iyr]{S)Sd-^ > 0. 

Remark 2.1 Since p has the meaning of the difference between fluid and external pressures, it can 



be negative. In this case, more than one value of u may satisfy (12). To ensure that the coupled 



model is not defective, one has to prescribe a particular rule for choosing the root of the cubic 
equation ( 12 1. In our numerical experiments, we take u which is the closest to |ro„j|~^ J u • nds. 



Remark 2.2 The boundary condition ( 12 ) is the combination of fluid and energy fluxes and so it 



does not guarantee to conserve the 'mass' of the entire coupled system. Although, we do not observe 
any perceptible generation or loss of mass in our numerical experiments, it does not necessarily 
mean that for all problems this effect should be negligible. Actually, one may consider any other 
linear combination of fluid and energy fluxes coupling on 3D-1D boundary to compromise between 
energy stability and mass conservation. 
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Figure 2: The schematic coupUng of ri"^, f^sD and JI^'d"" domains 



In practice, one may also be interested in coupling the ID fluid model to the upstream boundary 
of the 3D domain. Hence, we now consider three domains f^"}^, ^3d, and fi^'o™! as shown in 
Figure [2[ In il^^ and fi^g™" the simplified ID model is posed and in the full three-dimensional 
Navier-Stokes equations are solved. The domain fi"}?, is coupled to the inflow (upstream) boundary 
of r^sD and JI^'eT" coupled to the outflow (downstream) boundary of fiaD- The downstream 
coupling is described above. In the literature, it is common not to distinguish between upstream and 
downstream coupling boundary conditions. For example, in [51 133j conditions (|8|,([9]) are assumed 
both on upstream and downstream boundaries. Following this paradigm, one may consider ([8]),(|9]) 
or energy consistent conditions Q, ( 12 1 as the coupling conditions on Fin between ID model posed 



in r2-^p and 3D model posed in ilsD- Note that in entirely 3D fluid flows simulations, inflow and 
outflow boundary conditions usually differ. If a numerical approach to 1D-3D problem is based on 
subdomains splitting (see the next section for an example), then it is appropriate to distinguish 
between upstream and downstream coupling conditions. Thus, we impose the upstream coupling 
conditions in such a way that the 3D problem is supplied with the Dirichlet inflow boundary 
conditions. This is a standard choice for incompressible viscous fluid flows solvers and is especially 
convenient if third parties or legacy codes are separately used to compute 3D and ID solutions, 
and they communicate only through coupling conditions. 

For the upstream boundary, Fin , we introduce a reference velocity profile g(x), x e Fin, such 
that Jp g • nds = 1. Then the boundary condition on Fin is Dirichlet, given by 



Setting a 
equation 



Uin -ag on Fin. (14) 
Su\x=b ensures the continuity of the flux (|9| on Fi„. If a is found to satisfy the 



p{b)a + ^a^ J |g|2(g • n)ds = {pSu + ^Su'')U^t 



then the coupling condition ( 12 ) is valid on Fi„. Two more scalar boundary conditions are required 



for the ID model in r^"^. We assume that u or p are given in x = a and in a; = an absorbing 
condition is prescribed: in computations we set {Su)x = in a; = 6; another reasonable absorbing 
condition would be setting the incoming characteristic equals zero. On the downstream end of 
^ieT"; we also assume an absorbing boundary condition. 

We summarize the properties for the 3D-1D coupling introduced in this section: 



• It ensures the energy balance, as stated in Theorem |2.1[ 

• The inequality ([2| is not assumed; 

• It can be easy decoupled with splitting methods into the separate ID problems and the 3D 
problem with usual inflow-outflow boundary conditions on every time step. 
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3 Discretization and algebraic solver 



In this section, we introduce a splitting numerical time-integration algorithm based on sub- 
domain splitting. Further, we consider a fully discrete problem and review one state-of-the-art 
algebraic solver. 

Denote by 72",^", S*", u", and p" approximations to the corresponding unknown variables at 
time t = tn. Given these approximations, we compute u""^^ , p"''^^ , S"~^^, u"+^, and for 
t = tn+i {At — tn+i — tn) in three steps: 

Step 1. Integrate ([t]) for t g [i„,i„-|-i], with given u{tn+i) on the upstream end of the interval 
17"^ and the absorbing downstream condition at x = 6. 



Step 2. Set Ui„ according to ([Ti]), with {u,p,S} = S-^+i}, and compute p* and S* 

as the linear extrapolations oi p\x=d a-nd S\x=d from times i„ and tn-i- Solve the linearized 



Navier-Stokes problem in fisD for u"+^ ^n+i. 



1 

2At 



(3u 



n+l 



4u" + u"~^) + (2u" - u"-^) • Vu 



ri+l 



n+l 



Vp 



n+l 



cn+l 



( 



9u 



n+l 



dn 



div u 



P"+'n||r„ 



n+l _ 



= 0, 



p n. 



Step 3. Find u^+%=d from 



u 



n+l 



• nds 



u 



n+l 12/' n+l 



(u"+^ • n)ds. 



Now, using m""*"^ for boundary condition in x = d and the absorbing boundary condition in 
a: e, we integrate for t G [i„,i„+i] to find S'"+i in 1^5'°™. 

For the numerical integration of the ID model equations, we use a first order monotone finite 
difference scheme applied to the characteristic form ([T]), see [35!. To handle the 3D model, one 
has to solve on every time step the linearized Navier-Stokes equations, also known as the Oseen 
problem: 

- i/Au -t- (w • V) u + Vp = f 

in risD, 

divu = (15) 
u|r,„uro = g, ~ pn)|r„„, = 

where u — u^^^, p = p""*"^ — p*{d), the body forces term and the advection velocity field depend 
on previous time velocity approximations, f = (2At)~^(4u" — u"~^), w — (2u" — u"""'^), and 
P = 3(2Ai)~^. Here and in the remainder of this section, we dropped out the time-step dependence 
(n+l) index for unknown velocity and pressure. 

To discretize the Oseen problem (15), we consider a conforming finite element method. Denote 
the finite element velocity and pressure spaces by V/j C i/^(r23D)'^ and Qh C L^(ri3D), respectively. 
Let be the subspace of Yh of all FE velocity functions vanishing at Fin n Fq. 

The finite element problem reads: Find Uh € Yh, u/i|ri„ = u^J^, and ph G Qh satisfying 



with 



a(uft,v;0 - (p,i,divv,i) + (g,,,divu,j) {ih,Vh) V v/^ e Y^, qh e Qh, 
a{uh, ^h) = P{uh, v,J -I- i^(Vuft, \/Vh) + (w • Vu,,., v,,). 



(16) 
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Let — for ip,(l) £ W^. We assume the ellipticity, the continuity, and the stabihty 

conditions: 

I3i\\vh\\^ < ah{vh,Vh), a(v?„u,i) < /32||vh||v||u,,||v Vv^^u^eV^, (17) 

2 II ||2 ^ {qh, divvhf ^ ^ na\ 
7i MhW < sup V qh e Qh, (18) 

{qh, div V h) <-/2\\qh\\\Wh\W V G Q,,, v,, G V^, (19) 



with positive mesh- independent constants /?2, 7i, and 72. Condition (18) is weU- known as the 
LBB or inf-sup stabihty condition [9]. 

Let {4>i}i<i<n and {ipj}i<j<m be bases of and Qh, respectively. Define the fohowing 
matrices: 

Aij = a{(f)j,(f>i), Bi^j = -(div(/)j,'0i)- 



The hnear algebraic system corresponding to (16) {the discrete Oseen system) takes the form 

The right hand side {f,g)^ accounts for body forces and inhomogeneous velocity boundary con- 
ditions. To solve ( pO| ), we consider a Krylov subspace iterative method, with the block triangular 
preconditioner |12[ I21| : 



(21) 



The matrix ^ is a preconditioner for the matrix A, such that A ^ may be considered as an inexact 
solver for linear systems involving A. The matrix S* is a preconditioner for the pressure Schur 



complement of (20 1, S = BA ^B^ . In the algorithm, one needs the actions of A ^ and S ^ on 
subvectors, rather than the matrices A, S explicitly. Once good preconditioners for A and S are 
available, a preconditioned Krylov subspace method, such as GMRES or BiCGstab, is the efficient 
solver. In the literature, one can find geometric or algebraic multigrid (see, e.g., |12j and references 
therein) or domain decomposition p7l l3T] algorithms which provide effective preconditioners A for 
a range of v and various meshes. We use one V-cycle of the algebraic multigrid method [30] to 
define A"^ . 

Defining an appropriate pressure Schur complement preconditioner 5*^^ is more challenging. 
In this paper, we follow the approach of Kay ct al. f5(T. First, we define the pressure mass and 
velocity mass matrices: 

{Mp)ij = {ipj,ipi), {Mu)ij = {(j)],(l)i). 

The original pressure convection-diffusion (PCD) preconditioner, proposed in [20], is defined through 
its inverse: 

:= M-^ApL-\ (22) 

Here denotes an approximate solve with the pressure mass matrix. Matrices Ap and Lp are 

approximations to convection-diffusion and Laplacian operators in Qh, respectively. Both Ap and 
Lp (explicitly or implicitly) assume some pressure boundary conditions to be prescribed. 

If Qh defines continuous pressure approximations, one can use the conforming discretization of 
the pressure Poisson problem with Neumann boundary conditions: 

{Lp),,j = (V7/'j,V^,)- 

Likewise, Neumann boundary conditions are conventionally used to define the pressure convection- 
diffusion problem on Qh- However, the optimal boundary conditions setup both for Lp and Ap 
depends on the type of the boundary and flow regime, see [TJ] [5S| . 
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We use a modified PCD preconditioner defined below. Tliis modification partially obviates 
the issue of setting pressure boundary conditions and is consistent with the Cahouet-Chabard 
preconditioner [TU], if the inertia terms are neglected. The Cahouet-Chabard preconditioner is the 
standard choice for the time-dependent Stokes problem and enjoys the solid mathematical analysis 
in this case [24]. To define the preconditioner, we introduce the discrete advection matrix for 
continuous pressure approximations as 



Then the modified pressure convection-diffusion preconditioner (mPCD) is (compare to (22) ): 



S-i ^M-^ + {131 + Np){BM-^B^)-\ 

where is a diagonal approximation to the velocity mass matrix. 

Regarding the numerical analysis of the algebraic solver used here, we note the following. The 
eigenvalues bounds of the preconditioned Schur complement: 

< ci < |A(S'S^-i)| < Ci, (23) 

were proved for /3 = and the LBB stable finite elements in ^13) and for a more general case 
in |25j . The constants ci, Ci are independent of the meshsize h, but may depend on the ellipticity, 



continuity and stability constants in (17|-(19), and thus may depend on the problem parameters. 
In particular, the pressure stability constant 71, and so ci from ( [23^ , depends on the geometry 
of the domain D, [11 (tending to zero for long or narrow domains) and for certain FE pairs 71 
depends on the anisotropy ratio of a triangulation [T. Both of this dependencies require certain 
care in using the approach for computing flows in 3D elongated domains with thin and anisotropic 
inclusions (prototypical for simulating a flow over IVC filter). 

Characterizing the rate of convergence of nonsymmetric preconditioned iterations is a difficult 
task. In particular, eigenvalue information alone may not be sufficient to give meaningful estimates 
of the convergence rate of a method like preconditioned GMRES [18]. Nevertheless, experience 
shows that for many linear systems arising in practice, a well-clustered spectrum (away from zero) 
usually results in rapid convergence of the preconditioned iteration. This said, we should mention 



that a rigorous proof of the GMRES convergence applied to ( 20 ) , with block-triangular precondi- 



tioner (21 1, is not available in the literature (except the special case, when 5* is symmetric [?T1I5]). 



Thus, the numerical assessment of the approach is of practical interest. 

4 Accuracy and efficiency assessment 

In this section, we validate the accuracy and stability of the solver for the 3D- ID coupled fiuid model 
by (i) comparing the computed discrete solutions against an analytical solution for a problem with 
simple geometry; (ii) computing the drag coefficient and the pressure drop value for the flow around 
the 3D circular cylinder. The Taylor-Hood P2-P1 elements were used for the velocity-pressure 



approximation. The resulting linear algebraic systems (20 1 were solved by the preconditioned 
BiCGstab method. The initial guess in the BiCGstab method was zero on the first time step and 
equal to the (u",p") for the subsequent time steps. The stopping criteria was the 10~® decrease 
of the Euclidean norm of the residual. 

4.1 Test with analytical solution 

First we consider an example with analytical solution. The 3D domain is fl^i^ = {x G M'^ | x £ 
(— 1, 1), y"^ + < 1}. The circular cross sections are the inflow and outflow boundaries. Domains 
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max 1 u — u/J r2 
te[o,T] 


(/„^||V(u-^,)|pdt)^ 


(^^||^.-«,|pdt)^ 




mesh 1 


0.37 


0.30 


0.30 


8.10E-2 


mesh 2 


8.60E-2 (2.11) 


6.13E-2 (2.29) 


5.97E-2 (2.33) 


4.04E-2 (1.00) 


mesh 3 


2.66E-2 (1.69) 


1.68E-2 (1.87) 


1.62E-2 (1.88) 


2.02E-2 (1.00) 



Table 1: Errors to analytical solution in fisD on the sequence of refined meshes. Reduction orders 
are given in brackets 





in 


"id 


in nfg"" 




(^^||^.-u,fdt)4 


ij^ \\s ~ s^rdt)'2 


(/„^iiu-u,fdt)4 {j^\\p-p,rdt)'^ 


mesh 1 
mesh 2 
mesh 3 


0.25 

5.78E-2 (2.11) 
1.43E-2 (2.02) 


8.84E-4 
2.03E-4 (2.12) 
5.01E-5 (2.02) 


0.19 7.02E-4 
3.67E-2 (2.37) l.OOE-4 (2.81) 
5.64E-3 (2.70) 2.95E-5 (1.76) 



Table 2: Errors to analytical solution in il^^ and il^j'^*" on the sequence of refined meshes. 
Reduction orders are given in brackets 



fJ"^ and ri^*™" are two intervals of length 5. The analytical solution is given by 

5 = cos(27ri) + S^-l, u = l-cos(27rt), p^c^f{S), in J^^'g U 17^^™, 

2S \^ (24) 

— (l-cos(2^t))(l-y2-z2),0,0j , p=10{l-x)+p, in rigD, 

with S — TT, p ~ 1, c — 350, v — 1. This solution satisfies the continuity of flux condition on the 
coupling boundaries. The right-hand sides ip, and f were set accordingly. In this test, the 3D 
domain was triangulated using the global reflnement of an initial mesh, resulting in the sequence 
of meshes (further denoted by mesh 1, mesh 2, mesh 3), with the number of tetrahedra Ntet = 
1272, 8403, 63384, respectively. Since we use the first order scheme for the ID problem, the mesh 
size in ^l^^ and fi^g"" was divided by 4 on each level of refinement: Ax = 5/16, 5/64, 5/256. The 
corresponding time step was halved for every spacial refinement, so we use At = 0.02,0.01,0.005 
for mesh 1, mesh 2, and mesh 3, respectively. 

Based on the energy balance ( 13 1, the natural norms for measuring error in are C(0, T, L^{il3D)) 



and L^{0, T, H^{fl3u)) for velocity and L^{0, T, i'^(5l3D)) for pressure. These norms and, addition- 
ally, i^(0, T, L^{fl3Y))) for velocity error are shown in Table [ij The error norms in the ID domains 
coupled to the inflow and the outflow boundaries of ft^D are shown in Table [2] We observe the 
expected second order of convergence for all variables except for pressure in 3D. We remark that 
the integral in time error norms were computed approximately using the quadrature rule: 

/ ||V(^.-u,)||i.dt«At^||V(^.(nAt)~<')|li.(^3^), N^T{At)-\ 

Other integral norms were computed in the same way. 

The average number of iterations of the BiCGstab method for solving ( 20 ) , with block-triangular 
preconditioner (21 1, are shown in Tables [3] and |4] Table [3] shows that the convergence of the 



preconditioned BiCGstab method depends only slightly on the viscosity parameter and improves 
when the grid is refined. The latter observation is consistent with /i-independent eigenvalue bounds 



in ( 23 1 and with numerical results reported in jl2j for steady problems. Such robust behavior with 
respect to is observed only for sufficiently small values of the time step At. The results in Table |4] 
show that for small v-s the preconditioned BiCGstab method fails to converge unless At = 0.01. 
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v = l 


V = 0.1 


ly = 0.01 


ly = 0.001 


mesh 1 


10.7 


10.7 


13.4 


14.5 


mesh 2 


5.59 


6.69 


8.09 


8.99 


mesh 3 


4.42 


4.42 


6.50 


7.88 



Table 3: Average number of iterations of the preconditioned BiCGstab iterations for different 
meshes and viscosity parameters. 



At 


1 


0.5 


0.1 


0.05 


0.01 


0.005 


0.001 


0.1 


11.78 


10.78 


14.00 


* 


* 


* 




0.05 


7.16 


6.84 


7.21 


12.00 


* 


* 




0.01 


4.28 


4.45 


5.65 


6.27 


6.30 


6.63 


7.14 



Table 4: Average number of iterations of the preconditioned BiCGstab iterations for varying time 
step and viscosity parameters. The results are shown for mesh 2. 

4.2 Flow around circular cylinder 

Motivated by the simulation of blood flows over an intravenous filter, we experiment with flows in 
a 3D domain having an inclusion and coupled with ID model at the outflow boundary. Interesting 
statistics for such applications are the drag force acting on inclusions and the pressure drop. To 
validate the ability of the 3D solver to predict these statistics, we consider two benchmark problems 
of channel flows past a 3D cylinder with circular cross sections P^l, ^ . The 3D flow domain with 
a cylinder is shown in Figure |3] The no-slip and no-penetration boundary condition u = is 
prescribed on the channel walls and the cylinder surface. The inlet velocity is given by 

Uin = {l6Uyz{H -y)iH~ z)/H\ 0, 0)^ on Ti^, 

here H = 0.41m is the width of the channel. The kinematic viscosity of the fluid in this test is 
V = lO^^m^/s and its density is p = Ikg/m?. The Reynolds number, Re = v^^DU , is defined 
based on the cylinder width D = 0.1m and ?7 = |J7. We consider the following two benchmark 
problems from . 29J : 

• Problem PI: Steady fiow with Re = 20 {U ^ 0.45m/s); 

• Problem P2: Unsteady fiow with varying Reynolds number for U — 2.25 sin(7rt/8)m/s, 
t e [0,8]. 

The benchmarks setups do not specify outflow boundary conditions. Hence, on the outflow bound- 



ary we apply the 3D-1D coupling using the new conditions (|8|), ( 12 ) so that numerical performance 
of the coupling can be verified. 

The statistics of interest are the following: 

• The difference Ap = p(x2)— p(xi) between the pressure values in points xi — {0.2, 0.205, 0.55} 
and X2 = {0.2,0.205,0.45}. 



• The drag coefficient given by an integral over the surface of the cylinder 5* C F 



wall ■ 



Here n = [nj-^ny^nz)^ is the normal vector to the cylinder surface pointing to VL and t 
(— Uz, 0, rij;)^ is a tangent vector. 
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Figure 3: The flow domain fiao in the benchmark problem of the channel flow past a 3D cylinder. 



Table 5: Problems PI: Computed and reference values of drag and pressure drop. 



mesh 


C*drag 


% orr 


Ap 


% err 


Niter 


coarse 


6.149 


0.58% 


0.1679 


1.81% 


11.5 


fine 


6.196 


0.17% 


0.1678 


1.87% 


10.5 


Schafer & Turek 


[6.05,6.25] 




[0.165,0.175] 






Braack & Richter 


6.185 




0.1710 







Table 6: Problems P2: Computed and reference values of drag and pressure drop. 



mesh 


^drag 


% err 


Ap(t = 8) 


Niter 


coarse 


3.273 


0.76% 


-0.115 


11.7 


fine 


3.311 


0.39% 


-0.107 


10.6 


Schafer & Turek 


[3.2,3.3] 




[-0.11, -0.09] 




Bayraktar et al. 


3.298 
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For problem P2, the reference velocity in ( [25| ) is U = 2.25m/s. 

For these benchmark problems, the papeFpS] collects several DNS results based on various finite 
element, finite volume discretizations of the Navier-Stokes equations and the Lattice Boltzmann 
method. In [55], the authors provided reference intervals, where the statistics are expected to 
converge. Using a higher order finite element method and locally refined adaptive meshes, more 
accurate reference values of Cdrag and Ap were found in |8] for the steady state solution (problem 
PI) and in [3] for unsteady problem P2. For the computations we use two meshes: a 'coarse' 
and a 'fine' ones, both adaptively refined towards cylinder. The coarser mesh is build of 35803 
tetrahedra, which results in 53061 velocity d.o.f. and 8767 pressure d.o.f. for the Taylor-Hood 
P2-P1 element. The finer mesh consists of 51634 tetrahedra, which results in 73635 velocity d.o.f. 
and 12321 pressure d.o.f. Both coarse and fine mesh consist of regular tetrahedra. The refinement 
ratio is about 20 and 60 for the coarse and the fine meshes, respectively. We remark that the fine 
mesh has four times as many tetrahedra touching the cylinder as the coarse mesh. The time steps 
are 5t = 0.002 and 6t = 0.001 for the coarse and the fine meshes, respectively. 



Problem P2: drag coefficient 




Problem P2: drag coefficient 



















Coarser grid 






Finer grid 

Reference 











4 

time 



Figure 4: Evolution of the drag coefficient for unsteady flow around cylinder: coarse and fine grid 
results and reference results. The right figure zooms the plot for time in [3.8,4.2]. 

We first show in Tables [5] and [6] results for problems PI and P2 obtained with the coarse and 
the fine meshes. For all settings, the computed values are within "reference intervals" from [29] 
(except C^^g for problem P2, but in this case the upper reference bound appears to be tough). 
The computed drag coefficients were well within 1% of reference values and pressure drop within 
2%. This is a good result for the number of the degrees of freedom involved. Indeed, the results 
shown in [SJ O [H] for meshes with about the same number of degrees of freedom show comparable 
or worse accuracy. In Figure [4| we show the computed evolution of the drag coefficient for problem 
P2 and compare it to the reference results. The computed drag coefficients match the reference 
curve very well. We conclude that the conforming finite element method with the coupling outfiow 
conditions is a reliable and stable approach for the simulation of such flow problems. 

5 Simulations of a flow over a model IVC fllter 

The development of endovascular devices is the challenging problem of cardiovascular medicine. 
One example is the design of vascular filters implanted in inferior vena cava (IVC) to prevent a 
blockage of the main artery of the lung or one of its branches by a substance that has traveled 
from elsewhere in the body through the bloodstream. The filter is typically made of thin rigid 
metal wires as illustrated in Figure[5](left). Numerical simulation is an important tool that helps in 
finding an optimal filter design. Thin and anisotropic construction of a IVC filter requires adaptive 
grid refinement and makes computations of flows in such domains not an easy task. In this section, 
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Figure 5: Left: An example of intravenous filter (Corned Co.); Right: ID inflow IVC waveform 
used in computations. It was designed by interpolating the IVC Dopplcr blood flow waveforms 
from [31]. 



we demonstrate the ability of the numerical method to treat such problems in a stable way. One 
statistic of interest here is the drag force experienced by a filter. We recall that in this paper we do 
not account for the elastic properties of the vessel walls, which are otherwise important in practice. 

We consider a segment (4.5cm long) of IVC with elliptic cross section 1.6 x 2.4cm . The filter 
is placed on the 0.5cm distance from inflow, it is 2cm long and the diameter of its 12 wire legs is 
0.5mm. Blood is assumed to be incompressible fluid with dynamic viscosity equal to 0.0055Pa x s 
and density equal to Ig/cw?. 

A blood flow in IVC is strongly influenced by the contraction of the heart. The IVC have 
pulsatile waveforms with two peaks and reverse flow 39 occurring on every cardiac cycle. We 
consider the Doppler blood flow waveforms of IVC reported in [3S] and approximate them by a 
smooth periodic function plotted in Figure [5] (right). Note that the presence of significant reverse 
flows in IVC differs this problem from computing arteria fiows, where such phenomenon does not 
typically occur. 

On the inflow and outflow, the 3D vessel is coupled to ID models as described in section |2.3[ 
Each ID model consists of equations ([3|-([5]) posed on intervals of 5cm length. Periodic velocity 
with waveform as shown in Figure [5] is prescribed on the upstream part of the ID model coupled 
to Fin. The maximum ID model velocity of 12cm/sec yields the maximum inlet velocity in fiso 
of about 24cm/sec. This agrees with the measurements in |39j . The coupling conditions are the 
same regardless of the mean flow direction. 

The mesh was adapted towards the fllter, so the ratio of largest and smallest element diameters 
was about l.le + 2, the maximum elements anisotropy ratio was about 14. The resulting mesh 
is illustrated in Figure [6j The time step in 3D model was set equal to O.OOlsec. The BiCGstab 
iterative method, with preconditioner (21 1 was used to solve discrete Oseen subproblems. The 
stopping criterion was the reduction of the residual by the factor of 10^. The average number of 
linear iterations on every time step was about 35. We found that choosing time step larger for 
this problem, leads to the significant increase of the linear iteration counts and makes 'long time' 
computations non-feasible. 

We visualize the computed solutions in Figure [7] by showing the values of the x-component of 
the velocity in several cutplanes orthogonal to x-axis. Behind the filter the velocity x-component 
eventually has negative values, indicating the occurrence of circulation zones and 'returning' fiows. 
Note that the solution behind the fllter is no longer axial-symmetric: a perturbation to solution 
induced by non-symmetric tetrahedral grid is sufficient for the von Karman type fiow instability 
to develop behind the filter. 
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Figure 6: The visualization of the adaptive mesh for the flow over a model IVC filter problem: the 
top-left picture shows the boundary surface triangulation; the top-right picture shows the cutaway 
views of the tetrahedral grid. The bottom picture shows the zoom of the mesh in the neighborhood 
of the filter's 'head'. 
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Figure 7: The visualization of the velocity x-component in several cutplanes orthogonal to x-axis 
for times t S {3.06s, 3.34s, 3.39s, 3.52s, 3.66s, 3.92s}. One may note the occurrence of 'returning' 
flows behind the filter even for 'forward' mean flow. 
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Figure 8: Left: The evolution ol the drag force {gsm/s^) for the JVC fifter. Right: The evolution 
of the mean axial velocity (sm/s) in the middle point of the ID model before and after the 3D 
domain with IVC filter. 

Figure [s] (left) shows the time evolution of the drag force experienced by the filter. After the 
instantaneous start, the flow needs few cycles to obtain the periodic regime. In general, the drag 
force follows the pattern of the inflow waveform. In particular, the filter experiences forces both 
in downstream and upstream directions at different periods of the cardiac cycle. The right plot 
in Figure |8] shows the mean axial velocity in the middle point of the ID model before and after 
the 3D domain with cava filter. It is remarkable that after few cycles, when the flow is periodic, 
the waveforms in the ID domains coupled to upstream and downstream boundaries are very close. 
This suggests that the coupling conditions are efficient in conserving averaged flow quantities such 
as mean flux. 

6 Conclusions 

We reviewed the 3D and ID models of fluid flows and some existing coupling conditions for these 
models. New coupling conditions were introduced and shown to ensure a suitable bound for the 
cumulative energy of the model. The conditions were found to perform stable in several numerical 
tests with analytical and benchmark solutions. For the example of the flow around IVC filter, the 
coupled numerical model was found to capture the periodic flow regime and correct ID waveforms 
before and after 3D domain. The model was able to handle 'opposite direction' flow, i.e. the 
flow where the 'upstream' boundary (boundary with Dirichlet boundary conditions) becomes the 
outflow boundary for a period of time. The preconditioned BiCGstab method with one state- 
of-the-art preconditioner applied to the linearized finite element Navier-Stokes problem performs 
well. However, often the time step should be taken small enough to make the linear solver converge 
sufficiently fast. Overall, the coupled 3D-1D model together with the conforming finite element 
method and preconditioned iterative strategy was demonstrated as a reliable tool for the simulation 
of such biological flows as the flow over an inferior vena cava filter. 
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